!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine model_variables (ro, mx, my, mz, bx, by, bz, en, gx, gy, po, x, y, nx, ny)
!==================================================================================================================================|

   use parameters
   implicit none

   integer(4), intent(in) :: nx, ny

   real(8), dimension(nx), intent(in) :: x
   real(8), dimension(ny), intent(in) :: y
   real(8), dimension(nx, ny), intent(inout) :: ro, mx, my, mz, bx, by, bz, en, gx, gy, po

   integer(4) :: i, j, j_total, zeroy

   real(8) :: current_sheet_center_x
   real(8) :: current_sheet_center_y
   real(8) :: current_sheet_len_x
   real(8) :: current_sheet_len_y
   real(8) :: pr, dx, dy, g0, te_chr, te_cor, w_tr, h_cor
   real(8), dimension(nx, ny) :: te
   real(8), dimension(ny_total_finest) :: te_total_y

!----------------------------------------------------------------------------------------------------------------------------------|
!  Chromospheric evaporation
!----------------------------------------------------------------------------------------------------------------------------------|
   dx = x(2) - x(1)
   dy = y(2) - y(1)

   g0 = -1.0
   te_chr = 1.0d0
   te_cor = 1.0d0
   h_cor = 10.0d0
   w_tr = 1.0d0

   current_sheet_center_x = 0.0d0
   current_sheet_center_y = 0.0d0
   current_sheet_len_x = 0.05 * rangex
   current_sheet_len_y = 0.05 * rangey

   if (.not. allocated (y_total)) allocate (y_total(ny_total_finest))
   if (.not. allocated (ro_total_y)) allocate (ro_total_y(ny_total_finest))

   if (ro_total_y_flag) then

      zeroy = (nint ((ny_total - 2 * ng) * zero_posy / rangey) + ng) * 2 ** (n_levels - 1) + 1
      dy = rangey / dble(ny_total - 2 * ng) / 2.0d0 ** (n_levels - 1)

      y_total(zeroy) = dy * 0.5d0

      do j = zeroy + 1, ny_total_finest
         y_total(j) = y_total(j - 1) + dy
      enddo
      do j = zeroy - 1, 1, -1
         y_total(j) = y_total(j + 1) - dy
      enddo

      do j = 1, ny_total_finest
         te_total_y(j) = te_chr + (te_cor - te_chr) / 2.0d0 * (tanh ((y_total(j) - h_cor) / w_tr) + 1.0d0 + eps)
      enddo

      ro_total_y(zeroy) = (2.0d0 * (te_total_y(zeroy - 1) + te_total_y(zeroy)) / 2.0d0 + dy * 0.5d0 * g0) /                        &
         (2.0d0 * te_total_y(zeroy) - dy * 0.5d0 * g0) * 1.0d3
      do j = zeroy + 1, ny_total_finest
         ro_total_y(j) = (2.0d0 * te_total_y(j - 1) + (y_total(j) - y_total(j - 1)) * g0) /                                        &
            (2.0d0 * te_total_y(j) - (y_total(j) - y_total(j - 1)) * g0) * ro_total_y(j - 1)
      enddo
      do j = zeroy - 1, 1, -1
         ro_total_y(j) = (2.0d0 * te_total_y(j + 1) - (y_total(j + 1) - y_total(j)) * g0) /                                        &
            (2.0d0 * te_total_y(j) + (y_total(j + 1) - y_total(j)) * g0) * ro_total_y(j + 1)
      enddo

      ro_total_y_flag = .false.

   endif

   do j = 1, ny
   do i = 1, nx

      te(i, j) = te_chr + (te_cor - te_chr) / 2.0d0 * (tanh ((y(j) - h_cor) / w_tr) + 1.0d0 + eps)

      do j_total = 1, ny_total_finest
         if (y_total(j_total) .ge. y(j)) exit
      enddo

      if (j_total .eq. 1) then
         ro(i, j) = ro_total_y(j_total)
      else
         ro(i, j) = (ro_total_y(j_total - 1) * (y_total(j_total) - y(j)) + ro_total_y(j_total) * (y(j) - y_total(j_total - 1))) /  &
            (y_total(j_total) - y_total(j_total - 1))
      endif
      if (abs (x(i) + 0.0d0) .le. 1.2d0 .and. y(j) .ge. 10.0d0) then
         te(i, j) = min (te(i, j) + t * 10, 50.0d0)
         ro(i, j) = ro(i, j) / 10.0d0
      endif
      pr = ro(i, j) * te(i, j)
      mx(i, j) = 0.0d0 * ro(i, j)
      my(i, j) = 0.0d0 * ro(i, j)
      mz(i, j) = 0.0d0 * ro(i, j)
      bx(i, j) = 0.0d0
      by(i, j) = 1.0d0
      bz(i, j) = 0.0d0
!      if (x(i) - current_sheet_center_x .lt. - current_sheet_len_x) then
!         by(i, j) = -1.0d0
!         bz(i, j) = 0.0d0
!      elseif (x(i) - current_sheet_center_x .gt. current_sheet_len_x) then
!         by(i, j) = 1.0d0
!         bz(i, j) = 0.0d0
!      else
!         by(i, j) = 1.0d0 * sin ((x(i) - current_sheet_center_x) / current_sheet_len_x * pi / 2.0d0)
!         bz(i, j) = 1.0d0 * cos ((x(i) - current_sheet_center_x) / current_sheet_len_x * pi / 2.0d0)
!      endif
      en(i, j) = (mx(i, j) * mx(i, j) + my(i, j) * my(i, j) + mz(i, j) * mz(i, j)) / ro(i, j) / 2.0d0 +                            &
         (bx(i, j) * bx(i, j) + by(i, j) * by(i, j) + bz(i, j) * bz(i, j)) / 2.0d0 + pr / gmm1
      gx(i, j) = 0.0d0
      gy(i, j) = g0
      po(i, j) = 0.0d0

   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
end subroutine model_variables
